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ABSTRACT 


The analysis and design of a submarine propulsor requires the ability to pre- 
dict the characteristics of both laminar and turbulent flows to a high degree 
of accuracy. This report presents results of certain benchmark computations 
based on an upwind , high-resolution, finite -differencing N avier-Siokts solver. 

The purpose of the computations is to evaluate the ability , the accuracy and 
the performance of the solver in the simulation of detailed features of viscous 
flows . Features of interest include flow separation and reattc.chment. sur- 
face pressure and skin friction distributions. Those features .uTt particularly 
relevant to the propulsor analysis. Test cases with a wiut range of Reynolds 
numbers arc selected; therefore , the effects of the connect* 1 ' t and the diffusive 
terms of the solver can be evaluated ^epautely. Test cases include flows over 
bluff bodies, such as circular cylinders and spheres, at various low Reynolds 
numbers , flows over q, flat plate with and without turbulence effects } and 
turbulent flows over arisymmeiric bodies with and without propulsor effects. 
Finally , to enhance the iterative solution procedure , a full approximation 
scheme V-cyclc multigrid method is implemented. Preliminary results indi- 
cate that the method significantly reduces the computational effort. 

ADMINISTRATIVE INFORMATION 

This investigation was authorized by the Office of .Naval Research under 6.2 Viscous 
Flow Program, in accordance with Program Element 62323N, Task Area R2332MS3, 
and Work Request Number N0O1495WX2OO47/AJ. This work was performed at the 
Naval Surface Warfare Center, Cardevock Division (NSWCCD), David Taylor Model 
Ba^in (I)TMB) under Work Unit 5060-567. 

INTRODUCTION 

Hydrodynamically speaking, propulsor components are lifting bodies that provide 
thrust for propulsion. To improve the propulsive performance, it is desirable to have; a 
lifting body with an optimum lift to drag ratio. Successful analysis and design requires 
the ability to predict the hydrodynamic forces on the lifting body, such as lift and drag, 
to a high degree of accuracy. Recently, large research efforts on computational fluid 
dynamics (CFD) have been directed to achieve such goals. The most practical approach 
is to derive some numerical methods for the Reynolds- aver ged Navier-Stokes equations 
(R.ANS). The more popular methods, depending on the manner the convective lerms 
are formulated, include the central differencing and the upwind differencing schemes. 
Both upwind (DTNS Code 1 ) and central (1FLOW Code 2 ) differencing formulations have 
been used for hydrodynamics simulations at the David Taylor Model Basiri (DTMB). 
As part of DTMB propulsor research efforts, a variation of the upwind schemes has 
been employed to simulate the flow through the blade rows of a turbomachinery 3, 4 . 
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For high Reynolds number flows, the Navicr-Stokes equations become convectively 
dominated and are hyperbolic in nature. Under such circumstances, the upwind differ- 
encing approach becomes attractive and has certain advantages. For a one* dimensional 
case, it can be shown that the information at any point propagates in the direction 
according to the sign of the local eigenvalues of the flux Jacobian. Consequently, a 
numerical procedure with the upwind differencing technique based on the direction of 
local wave propagation is physically more adaptive to the characteristics of the equa- 
tions. Besides, in the upwind scheme, the dispersive and dissipative errors are more 
closely balanced than in other schemes of equal or even highei order that use the same 
set of nodal points, regardless of the direction of the convection 5 . Tbe flux-differencing 
splitting approach suggested by Roe 6 is a popular upwind scheme for solving the in- 
compressible Navier -Stokes equations, since it docs not require the inviscicl flux to be a 
homogenous function of order one. Based on Roe’s approach 6 the flux-difference is split 
according to an approximate solution of the Riemunn problem. The solution provides 
the information about the direction of the wave propagation, which, in turn, is incor- 
porated into the discretized system to form an upwind scheme. The accuracy of the 
solution can be promoted to a higher order by reconstructing the primitive variables or 
the fluxes midway between two nodes with an extrapolation technique. In conjunction 
with the reconstruction process, the slope limiter or the flux limiter can be implemented 
to obtain total variation diminishing (TVD) property 7 . Viscous terms can be discretized 
with a second order central differencing scheme. A numerical formulation, based on the 
principles described above, aJlows the discontinuity of the solution to be resolved ove. 
only two adjacent nodes without causing the non-phv&ical oscillations. Therefore, this 
approach is sometimes called high resolution. 

This report presents the results of a numerical study designed to evaluate the ability 
arid accuracy of ctn upwind scheme in predicting certain flow features that are relevant 
to propulsor analysis. Some of the features of interest are flow' separation and reattach- 
ment. surface pressure and skin friction distribution. Test cases for present numerical 
study are selected so that (1) a wide range of Reynolds numbers are covered, (2) the 
boundaries are smooth and the distortions; of the grid systems are minimized, (3) the 
flow features are relevant to the proplusor analysis and (4) the test data are well ana- 
lyzed. Test eases include flows over bluff bodies, such as circular cylinders and spheres, 
at various low Reynolds numbers, flow's over a flat plate with and without turbulence 
effects, and turbulent flows over axisymmetric bodies; with and without propulsor ef- 
fects. To enhance the iterative solution procedure, a full approximation scheme (FAS) 
V- cycle multigrid method is implemented. A fast convergence rate is achieved as a 
result. 
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THE GOVERNING EQUATIONS 


The three-dimensional incompressible RANS equations based on primitive variables 
are formulated iu a boundary-fitted curvilinear coordinate; system. Using Chorin's arti- 
ficial compressibility formulation, 8 the* incompressible Navier-Stokes equation is written 
in conservation form for three-dimensional flow in Cartesian system as 

Qt + (£• - E;) x -f (F- - f;) y + (<?* - Gi) z 0 . (1) 

In Eq. 1, the dependent variable vectoT Q is defined as Q —■ (p, v , v , w) T , and the inviscid 
flux vectors E",F’ i and G\ the viscous shear flux vectors i'/'.-F*, and G’ v are given by 

E * “ (/?u.u J 4p, uv.uui) 7 

F~ = (Pv,uv,v 2 + p.v w) T 
G" —■ (8u;,uw, vw.w 2 4 p) 7 
E * =: J Re'' 1 (0,2u x ,u ! , 4 r x ,u- 4 w x ) T 
F’. =• Re~ } ( 0, u y -1- v x . 2v y , v 2 + !!;/ 

6" =•• Re" x (0, u 2 4 w x , u* 4 2u;- ) T . 

The coordinates x,y,z are scaled with an appropriate characteristic length scale L. The 
Cartesian velocity components u.v, w are nondimensionalized with respect to the free 
stream velocity, K*.. The normalized pressure is defined as p = (p - Po&)/pV^. The 
kinematic viscosity, r\ is assumed to be constant, and the Reynolds number is defined as 
Re = V x L/v. The artificial compressibility parameter, /?, monitors the error associated 
with the addition of the unst eady pressure term dp/dt in the continuity equation. The 
unsteady pressure term is needed for coupling the mass and momentum equations to 
make the system hyperbolic. 

Equation 1 can be transferred to a. curvilinear, body-fitted coordinate system (( ,£ , 
r/ ) through a coordinate transformation of the form 

( -■ ((x,y,z), < = £(x,y,z), and rj = i]{x . y, z). 

Thus, Eq. 1 becomes. 


(Q/J) t 4 (E - E v ) c 4 (F - F v ) ( 4 (G - G„) n = 0 


(2) 


with 

(E,F.G) t = { \T) [E’/J ,F * / J ,G'JJ) T ) 

and 

(E Ul F v ,G\) r = [T] (E'JJ , F;/J ,C:./J) t 
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where 


m= 


Cr Cy Cx 
£r £y £z 


L >/* >h V* J 


and the Jacobian of the coordinate transformation is given by 


./ 


-i 


det 


x < V< z ( 
*(. Vi H 

*i» Vv . 


The Jacobian is the ratio of volume elements in the two coordinate systems. For a proper 
transformation, neither J nor its reciprocal is zero. At present, the transformation is 
chosen so that the grid spacing in the computational domain is uniform and of one unit 
in length in all three spatial directions. The Cartesian derivatives of the shear fluxes are 
obtained by expanding them using chain rule expansions in the and rj directions. 


NUMERICAL DISCRETIZATION 

Let b j,k denote the integer indices of a grid point in the curvilinear system 
whose Cartesian coordinates are x, y,z. Each grid point serves also as the centroid 
of a control volume whose six bounding surfaces are formed by bisecting the distances 
between the centroid and its six adjacent grid points. The Cartesian flow variables such 
as u y v s w arid p are placed at each, grid point. The indices of the grid point are used 
as subscripts for the variables to indicate the association. To avoid introducing any free 
stream error, the metric terms such as and etc., are computed from r,t/,z 

data by using a second-order central-diffcrencing approximation of x ^ and x Ti etc., 

as described by Pulliam and Stegerh 

INVISC1D FLUXES 

At present, discretization of the inviscid fluxes of Eq. 2 is achieved by applying 
the Riemann solver to each direction of the coordinate system as suggested by Yee, 
Warming and Harten 10 . Consider a one-dimensional hyperbolic system of conservative 
laws 


(j),+ . ¥.«?)=», = ® 

where & — (, £ or ?/, and D(Q) is the Jacobian Matrix. Both Q and H{Q) arc column 
vectors with four components and D(Q) is a four by four matrix. 

Let the right eigenvectors be the column elements of matrix R and eigenvalues 
Aj, A 2 . A 3 and A 4 be the diagonal dements of matrix A, then the relationship between 
D,i? and A is given by a similarity transformation D - fi’Af? -1 . The row elements 
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of the matrix R~ } give an orthonormal set of the left eigenvectors. In the discretized 
system, the state Q/ at grid point / is considered as an averaged value in an interval, 
that is 

Q, = Ar 1 / 7 Q do. (4) 

Roe’s flux difference splitting is constructed by forming a mean value Jacobian D(Qu 1. Qt) 
such that 

H(Qt 4-1) - H(Q,) = D(Q w ,Q t ) (Q l+l - Q t ), 

and 


D(Qi+i-.Qi) - D(Qu. u Q,) as Qi+i -•» Qi . 

to evaluate the mean value (the local frozen value ) Jacobian at the interface ! ± J, the 
simple average value of the state Q at two adjacent grid points are used, that is 

A±i ~ (5) 

where 

Qi±t = ${Qt±Qi*i) 

By the similarity transformation, D l+ 1 can be written as 
•^t+j ~ (RRR-i )f+| 

” D l+l ~ D u{ ’ 

where A 4- and A - ' are the absolute values of the positive (speed of the right travelling 
wave) and the negative (speed oi: the left travelling wave) eigenvalues, respectively. 
With the mean value Jacobian locally defined* Eq. 3 can be decoupled into a system of 
four scalar equations with the eigenvalues representing the wave speeds of the Riemann 
problem. Let’s define the local characteristic variables W as the projection of Q into 
the left eigenvectors Z?** 1 ; therefore , W = R~ X Q. For any given two states Q\ and 
Qi+i, the flux at the interface can be expressed, in term of flux difference, as either 


H l+ i *= H, + RA-A l+ iW 

(6) 

H l+ ^Hu i - RA + A (+ .1V , 

(V) 
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where A l+ iW = (W, +l - W,)/2 . 

Assuming that H, R and A are constant, with Eqs 6 and 7, Eq. 3 can be written as 
a system of four scalar hyperbolic conservative equations lor the characteristic variables 
IF , that is, 


/ H 'm \ 

J f ~ - 0 for 7?? = 1,2. 3,4. 


(S) 


An implicit delta form of Eq. 8 is 


1 


\~J + (A m )j_lA;_.l 

= (A m )/ + i A j+ i 14 7 t " - (A 4 ),_i A ; _i W™ , 


AH- 


(9) 


where r is the time step size, n is the time step number, and AW- = If- 41 — W”. 

When is sufficiently differentiable, the local jump in W” at the interfaces l±K 
known as 




are replaced wdth 



( 10 ) 


(H) 


Equation 10 is only a first-order one-sided approximation of Eq. 11. To enhance the 
accuracy, the following relationships from a scalar scheme can be extended to a constant 
coefficient system by applying them scalarlv to each of the m scalar characteristic 
components in Eq. 8, 


+ 


and 


* 2 

H<> + + a - <*xnj - *t j.. )a,.i h" 

= K- - "=,) 

[<’ +")(4 , ;.J - »J;.I-|)A,. S U" + (1 - 


.(12) 


for 7n=l,2,3.4, with 
= , 


■1*3; 


(M) 
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and 


(A H W"/ A J+i W2)* for A ;± 1 IV” # o 
0 for = 0 


The order of the accuracy in the spatial derivative's presented above is determ ineci by the 
values of u). For ~ —1, the scheme is fully upwind second-order accurate. For \ * 

the scheme is upwind biased third-order accurate. Function $ is called the limiter; it 
is used to control unwanted oscillations in numerical schemes. Various designs of the 
limiter were found and successfully tested. 3 - 13 

Substitute Eqs. 12 and 13 into Eq. 9, and multiply Eq. 9 by the set of right eigen- 
vectors R from the left. A conservative high-resolution scheme for the nonlinear system 
is derived: 






A Q n = 


(RK-' R.-'K +i A ul Q* - nr >)JLi A,_xQ* 


(16) 


with 



and 


The right hand side of the Eq. 16 is evaluated at time level n; it is the spatial 
derivative of Eq. 2. and is designated as residual. Equation 16 represents the relationship 
between the residual at n ih time step and the correction of the solutions at n 4 ] time 
step). The correction and the residual approach to zero as the solutions approach to 
their steady state values. 

VISCOUS FLUXES 

The viscous fluxes in Eq. 2 are evaiuted by a second-order central differencing 
scheme. The computation of the fluxes require all nine metric coefficients at each 
of the six bounding surfaces of each computational cell. 


SOLUTION ALGORITHM 

Equation 16 can be extended to three-dimensional applications with the operator 
split method. The differencing schemes described previously are applied to each co- 
ordinate direction (. £, and rj independently; a summation over all directions gives 
the discretized approximation of a multi-dimensional flow problem. Upon forming the 
Jacobian matrices A. B, and C from invicid fluxes E, F, and G and A', Y . and 
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Z from the viscous fluxes F v . and G v , a finite-differencing form of Navier-Stokes 
equations can then be written as: 

{7j ~ + + 

- <B- + V)* j A*, +(£♦■+ »-),_, Vj 

- (C-+/), tl a 1+i +(C- + Z) t .jA,.j)*A<3*--fifiS(<?») . (17) 

Equation 17 is solved by an implicit hybrid algorithm where a symmetric planar 
Gauss;- Seidel relaxation is used in the ( direction and approximation factorization is 
applied in the £ and r/ directions: 


M - (ir + F),. +i A, 4 : + (S + + Y) H * H ' 
[M - (C- + Z) k+ iA k+ i + (C + + Z) k -iA k _i 


Q u+l = Q» + A Q n , 


AQ = -RES(Q”. C?" 41 ) , (18) 

A<? = MAQ , (19) 

(20) 


with M — I !(tJ) + (A~ -f A')- + i + (A 4 + where AfiS((2 r, % Q T ‘ 41 ) indicates the 

nonlinear updating of the residual while sweeping in the ( direction. 


MULTIGRID METHOD 

For certain simulations, in order to obtain meaningful results, a large number of 
finely sized grid points is; needed. The adverse effect of such a grid system upon the 
computation effort is that the rate of convergence deteriorates significantly. The analysis 
provided by Brandt 12 suggests that the low-frequency components erf the errors are 
efficiently approximated on coarse grids but are slow to convergence on fine grids. 
In addition, the high-frequency components must be approximated on fine grids. By 
utilizing interactively several scales of discretization, multigrid techniques resolve such 
conflicts and avoid stalling. 

To accommodate monlintearities, a full approximate scheme (FAS) is used. The 
discretized system of equations described previously can be represented as; 

!(<?) = -/*. ( 21 ) 

where L is the differencing, operator, Q is the unknown to be solved and i? is the 
residual- The iterative process will reduce the residual to zero as the steady-state 
solution is approached. The FAS procedures for solving F/q. 21 can be described as 
follows: 

(l) relax on the fine grid, Li k (Q^) = — R h , 
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(2) solve L 2 h(Q 7h ) 4 (/| A [L/>(<? , ‘)] - L 2 u(Jf l h Q h ) = -R? n on the coarse grid, and 

(3) replace Q h < — Q h + l 2>l (Q 7h — lJ^Q h ) on the fine grid. 

Tlie notation introduced here includes the difference operators at the fine grid Lf t 
and the coarse grid Lu . the restriction operators l\ h (for the approximation) and 
I% h (for the residual), and the interpolation operator 7^. 

RESULTS 

In the followings benchmark computations, in order to qualify the comparison be- 
tween the experimental measurements and computational resul ts, roof -mean-square 
(RMS) differences 2 are calculated. The RMS difference is defined as 

RMS= [l>\ c --<-)7A’]’ . 

where v r: and v m are computed and measured values respectively, N is the total 

number of data values used in the comparison, and the subscript i ranges from 1 to 
/V. Computations were carried out with 64-bit precision on a Silicon Graphics Power 
Challenge machine. 

TURBULENT FLOW ALONG A FLAT PLATE 

The structure of the; turbulent boundary layer along a flat plate has been investigated 
earlier by Ludweig and Tillmana 13 , and Wieghart and Tillmann 14 . It was found that 
the axial velocity profile in the inner onc-fifth of the boundary layer can be represented 
by the universal logarithmic law (excluding the viscous sublayer). The remaining outer 
four-fifths can be adequatedly expressed by the power law. Wieghart and Tillman n's 
data were collected in a wind tunnel test. The flow velocity was 33 rn/sec arid the 
average dynamic viscosity was 0.151 cm 2 /sec. Velocity measurements were taken at 
twenty -three locations ranging from 0.087 m to 4.987 in from the leading edge. The 
boundary thickness grew from 0.0335 cm to 0.9242 cm. The test data were compiled 
and presented at a 1968 Stanford turbulent flow conference 15 . 

In the present computation, the computational domain extends 8 m in stream wise 
direction. 0.16 m in cross flow direction, and 0.5 rn in the third direction: grid points 
used are 57., 61 and 5 in the respective directions. The Reynolds number (Re) is 2.2:<10 fi . 
Grid distribution in the cross-flow direction is non-uniform, and is clustered near the 
plate such that the t/ 4 coordinate of the first grid point off' the plate surface is less 
than 0.3. The Coiirani-Friedricbs-Lewy (CFL) number for present computations is 
10". For the turbulence modelling, the standard Baldwin and Lomax’s algebraic eddy 
viscosity formulation 10 is used. Figure la shows the skin friction coefficient C r along 
the* plate surface, and Figs, lb and 1 c show the axial velocity profiles at x=0.78 ni and 
4.98 m, respectively. The RMS differences indicate that the deviations between the 
measurements and the predictions are within the limits of the expected measurement 
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uncertainties. Figure 2 shows the velocity profile at x— 4.98 in presented in (y * , u + ) 
coordinates. It can bo seen that the predictions agree well with the measurer tents. In 
the turbulent zone where y* lies between 30 and 1000. both the predictions and the 
measurements fit the universal logarithmic law. However, none of the data lie within 
the sublayer. The convergence history presented in Fig. 3 shows that the residual value 
approaches the machine zero. 
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Fig. lb. Velocity profile at x=0.78 m. 


Fig. lc. Velocity profile at x=4.98 m. 


Fig. 1. Skin friction coefficient and velocity profile 
for turbulent flow along a flat plate. 
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Fig. 2. Velocity profile in (y + ,u + ) coordinates. Fig. 3. Conveigenee histories for computing 

turbulent flow along a flat plate. 


LAMINAR FLOW ALONG A FLAT PLATE 

For steady laminar flow along a flat plate with a- constant free stream velocity 
l‘o , the pressure gradient dp/ Ox in stream wise direction vanishes. The Navier-Stokes 
equations reduce to the Prandtl boundary-layer equations. A solution, known as the 
Bla-siiis solution after its originator, is obtained by assuming similar profiles along the 
plate .at every location along the plate. Blasius assumed that 



where y is the distance above the plate surface, ^ is the boundary layer thickness, and 
V V 

* *lJK ~ 9 ’ 

with x the distance from leading edge and R x the Reynolds number based on z. Under 
the similarity assumption, the Prandtl boundary layer equations can be further reduced 
to an ordinary differential equation. Solution can then be obtained numerically. In the 
present numerical simulation, the geometrical dimension used previously for turbulent 
flow is adopted and the Reynolds number is 3.64xl0 5 . The grid size is 129x129x3. 
figure 4a shows the skin friction coefficient along the plate surface, and Figs. 4b and 
4c show the profiles of axial velocity u/U 0 , and transverse velocity w/U a , respectively. 
The R.MS difference for each quantity is also presented. The results of the present 
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Fie. 4a. Skin friction coefficient. 


Similarity solution of axial velocity 

o Blasius solution 
Present computation 


5 qL Similarity olution of transverse velocity ; 


Blasius solution 


Present computation 


2-0 - RMS difference = 0.010 
in: 


570 0.2 0,4 0.6 0.8 1.0 

u/U 

Fig. 4b. Axial velocity profile. 


\jdr RMS difference * 0.023x10 
(w/UJx 10 3 

Fig. 4c. Transverse velocity profile. 


Fig. 4. Skin friction coefficient and velocity profiles 
for laminar flow along a flat plate. 


computation agree well with the solutions obtained by Biasings method. Part of the 
deviation can be attributed to the fact that the Blasius's solution is based on Prandtl 
boundary layer equations while the current computation is based on N&vier- Stokes 
equations. Figure 5 shows the convergence histories for solutions with multigrid (7 
ievels) and without multigrid (1 level) application. Considerable computing effort is 
saved with the application of the multigrid method. 
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Fig. 5. Convergence histories for computing 
laminar flow along a flat plate. 


FLOW PAST CIRCULAR CYLINDERS AT LOW SPEEDS 

For flow past a circular cylinder at Reynolds numbers (length scale is based on the 
diameter) below the critical value ( ~ 40) at which a Karman street is formed, the flow 
is steady and twin vortices exist behind the cylinder. There are many interesting flow 
related phenomena despite the simplicity of the geometry involved. The phenomena 
are the boundary layer development along curved surface, the flow separation, and 
wake reattachment. However, the details of these phenomena, such as the locations of 
the separation, the coordinates of the vortex center, and the wake’s length and shape, 
are Reynolds numbers dependent and pose great challenges for numerical simulation. 
For these and other reasons, this problem has drawn much attention in the past, both 
theoretically and experimentally. In light of the availability and quality of the data, 
this problem is selected as a benchmark case for present numerical study. 

The outer boundary of the computational domain h described by a circle with a 
radius 10 times that of the cylinder. An O-typc and orthogonal grid system is se- 
lected. The grid distribution in the radial direction is non-uniform and' is clustered 
near the surface. The distance between the cylinder surface and the first grid point is 
one- thousandth of the radius of the cylinder. For computation, symmetry is assumed 
and only the plane above the axis of symmetry is considered. The number of grid 
points in radial, circumferential and axial directions are 257, 257 and 3, respectively. 
The CFL number used for the following computations is 10 2 . The boundary conditions 
are: (1) a non-sup , non-penetrating, and vanishing normal pressure gradient on solid 
surface, (2) prescribed free stream values at; the upstream side of the outer boundary, 
(3) second-order extrapolations at the downstream side of the outer boundary and (4) 
periodicity in the spanwise direction. Jn order to describe the main features of the flow, 
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Fig. 6. Geometrical parameters of the closed wake. 


Coutanceau and Bollard's 17 geometrical parameters for the closed wake are adopted 
and shown in Fig. G. where parameters l and L are the width and length of the wake, 0 3 
is the separation angle and the vertex centers arc located on the (x>y) axes by a and 
b . It was observed by Taneda 18 that the twin vortices behind the cylinder appeared 
when the Reynolds number was greater than 6 and became unstable when the Reynolds 
number was greater than 45. Therefore, the Reynolds numbers 5 and 40 are selected as 
the lower and upper bounds for the present steady-state computations. Taneda*s 18 and 
Coutanceau and BouardV 7 flow visualizations were obtained by similar methods: by 
illuminating the light particles suspended uniformly in the liquid and by photographing 
in the direction normal to the lighted plane. Taneda 18 used aluminium and water in 
the tests, while Coutanceau and Bonard 37 used fine bright particles and the Vaseline 
oil ‘ MAR.COL 80 1 in their tests. The latter reference derived the particle velocity 
by measuring the length of the particle trajectory during the time of exposure. The 
reported inaccuracy was less than 2%. The wall influence was investigated by chang- 
ing the ratio A between the cylinder and the tank diameters 17 . Figure 7 shows the 
com puled wake shapes behind the cylinder at, different Reynolds numbers. Compared 
with the flow visualizations 17, lfi , the characteristics of the wake shapes near the sepa- 
ration and reattachment, points are well simulated. The result indicates that the twin 
vortices begin to develop at Reynolds number about 5 and it agrees with TancdaV 8 
observation. In Fig. 8. the currently computed separation angles at various Reynolds 
numbers arc compared with those computed earlier by Keller 19, 20 and Raal 21 arid those 
measured by Coutanceau and Bouard. 17 Figure 9 show's the relationship between the 
Reynolds numbers Ro and the wake lenglhs L(D , from both the current compulations 
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Reynolds number. Re Reynolds number, Re 

Fig. 8. Separation angles at different Reynolds Fig. 9. Wake lengths at different Reynolds 
numbers for flow past a circular cylinder. numbers for flow past a circular cylinder. 

and Tancda’s observations. 18 In Fig. 10, t he coordinates of the vortex center (a, b) (see 
Fig. 6) are plotted against the Reynolds number. Re. The experimental data were re- 
ported by Coutanceau and Bouard 17 with A , the ratio between the cylinder and the 
tank diameters, equal to 0.024. The small value of A indicates that the wall effect is 
relatively small. Figure 11 shows the similarity of the closed-wake shape. 17 For com- 
puted and measured wakes at different Reynolds numbers, when the wake width / and 
its distance from the rear stagnation point A' — /?. arc normalized with the maximum 
wake width / ma3: and wake length L> respectively, and then plotted against each other* 
the results merge into a single curve except at the regions near the cylinder wall. Fig- 
ure 12 show's the velocity similarity on the rear flow axis in the closed wake, where the 
velocity u is normalized with its maximum value u rnax . Figure 13 shows the compar- 
isons of the computed and measured velocities at rear flow axis at Reynolds numbers 
20 and 40. The RMS differences are comparable with the measurement uncertainty 
2%). Figure 14 shows the computed pressure distribution at the cylinder surface 
for Re=4(l- it is compared with the distributions observed at Rc=36 and Re— 45 by 
Thom. 22 and computed by Apelt 23 at Re=40. Thom’s approximate theory 24 for deter- 
mining the value of the pressure at the front stagnation point at low speed gives the 
result at Reynolds number 40, (1+7/Re) or 1.175, which agrees well with the value of 
the current, result 1.18. Figure 15 presents the convergence histories of the numerical 
simulation at Re=40. The 7-level multigrid solver improves the efficiency significantly. 
A fine grid size is needed for detailed computation, because with such a fine grid size 
the convergence is slow, The application of multigrid technique (7 levels) reduces the 
computational effort considerably for a given CPU time. 
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Fig. 12. Velocity similarity ori the rear flow axes 
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Fig. 13a. Re-" 20. Fig. 13b. Re=40. 

Fig. 13. Velocity distributions on the flow axis behind the circular cylinder. 
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FLOW PAST SPHERES AT LOW SPEEDS 

The flow and the solutions are assumed to be axisymmetric. The topology of the 
grid and the specification of the boundary conditions are the same as those used for 
coinouting the flow over the cylinder presented earlier* except that the reflective con- 
dition is applied in the circumferential direction. The solutions on different meridional 
planes aTe related by simple coordinate transformation. The CF'L number used for the 
computations is 10 2 . Figure 16 shows the computed wake shapes behind the sphere at 
various Reynolds numbers. The relationship between the Reynolds number Re and the 










I 


Fig. 16. Standing vortices downstream of a sphere. 
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wake length L is plotted in Fig. 17; also included are the test data from Tancda. 25 
The twin vortices behind the sphere appear at Reynolds number about 25. The surface 
pressure distribution at Reynolds number 100 is shown in Fig. IS; also included are the 
results from earlier computations. 26 2 ' The convergence histories for Re— 1 00 is shown 
in Fig. 19. Multigrid technique improves (he computational efficiency significantly. 



Reynolds number. Re 


Fig. 17. Wake lengths at different Reynolds 
numbers for flow past a sphere. 



Degree from front stagnation point Work unit 


Fig. 18. Surface pressure distribution F '£- 19. Convergence histories for computing 
for flow past a sphere at Re^lOO. flow past a sphere at Re=100. 
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TURBULENT FLOW PAST AXISYMMETRIC BODIES 

The purpose of this computation is to investigate the accuracy of the numerical 
scheme in predicting the characteristics of the flow in an expanding boundary layer 
under adverse pressure gradient at high Reynolds numbers. Axisymmctric bodies, des- 
ignated as DTMB-bodyl and DTMB-body2, were built and tested at the David Taylor 
Model Basin ''DTMBj 28 . The Reynolds number at the test condition was 6-fixlO 6 . The 
size of the C-type grid for the present computation is 197x3x146 in radial, circumfer- 
ential and axial directions, respectively. Modified Baldwin-Lom&x model 2 is used to 
simulate the turbulent flow. Figure 20 show's the computed (solid line) and measured 
(symbol) pressure distributions on surface of DTMB-bodyl. The RMS difference is 
0.019. The measurement uncertainty for pressure was ±0.015. It can be seen that 
the pressure gradient is zero at the parallel mid-body region, and adverse gradient is 
present at stern region. In Fig. 21 the computed velocity profiles at various axial loca- 
tions at stern region of DTMB-bodyl are compared with the measurntents. The RMS 
differences range from 0.012 to 0.025. The measure.!* ,t uncertainty for velocity was 
.1-0.025. In Fig. 22 the velocity profiles, at the parallel mid-body section, are plotted in 
the (y + ,u + ) coordinates. The value of the y+ of the grid nearest to the bodv surface 
is about 2.5. There are four grid points that lie within the laminar sublayer. The 
computed distribution of skin friction (solid line) on body suiface is shown in Fig. 23. 
Compared with the measurements, the RMS difference is 0.00029. The measurement 
uncertainty for skin friction was ±0.0002. The computed (solid line) and measured 
(symbol) turbulence shear stresses near the stern region at several axial locations are 
shown m Fig. 21; the RMS differences range from 0.008 to 0.015. The measurement 
uncertainty for turbulent shear stresses is ±0.01. 
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Fig. 20. Surface pressure distribution 
for DTMB-bodyl. 
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Fig. 21. Velocity profiles at stern region of DTMB-hodyl. 
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Fig. 22b. x/L=0.5765 
Fig. 22. Comparisons of wall law at mid-body section of DTMB-bodyl. 



Fig. 23. Skin friction coefficient 
along the surface of DTMB-bodyl. 
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Fig. 24a. x/L=0.755. 


Fig. 24b. x/L=0.S'J6. 




Fig. 24. Turbulent shear stresses near the stem region of DTMB-bodyl. 


24 



A propeller was placed on DTMB-bodyl whose center line is located at x/I.,=0.9S3. 
The ratio of propeller and body diameter is 0.54. Velocity measurements were taken 
at a distance of 0.227 propeller diameter upstream of the propeller. For numerical 
prediction, the propeller effect was simulated with a body force model. The computer! 

and measured velocity profiles under two different propeller operating- conditions are 
shown in Fig. 25. 

The same type of computations were carried out on DTMB-body2. The velocity- 
profiles at different axial locations are shown in Fig. 26. The RMS differences between 
the computed and measured values range from 0.008 to 0.058. 




Fig. 25b. Advance ratio=!.07. 


Fig. 25. Velocity profiles upstream of an operating propeller 
at different, advance ratios. 
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Fie. 26a. x/L=0.840. 



Fig. 26b. x/L=0.933. 
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Fig. 26. Velocity profiles at 
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CONCLUSIONS AND RECOMMENDATIONS 

An upwind, high-resolution, finite-differencing Navier-Stokes solver is used to sim- 
ulate steady-state incompressible flows with a wide range of Reynolds numbers. Ex- 
periments with measurements of high quality are selected as benchmark cases. The 
predictions are compared with the measurements by evaluating the root mean square 
(RMS) differences. In all cases, the RMS differences are compatible with the mea- 
surement uncertainties. For the low Reynolds number cases, the detailed features of 
the standing vortices behind the bluff bodies are successfully simulated. For the high 
Reynolds number cases, the skin friction coefficients, the structure of the turbulent 
velocity profiles and the turbulent shear stresses are correctly predicted. 

At low Reynolds Dumbers, the flows are dominated by diffusive process. The rate of 
convergence of fhe iterative procedure becomes: very slow', even for an implicit method 
with a rather high C'FL number. The situations can be improved significantly by 
implementing the multigrid technique. Compared with the single-grid approach, the 
multigrid solver is rather insensitive to the CFL numbers and an order of magnitude of 
Central Processor Unit (CPU) time is saved. 

In summary, the flow features that are relevant to the propulsor analysis, such as 
flow separation and reattachment, surface pressure and skin friction dstributions, can 
be predicted accurately and efficiently with an upwind RANS solver. The formulations 
of DINS 1 are similar to the formulations described in this report. It is expected that 
D’TNS 1 code may achieve the similar perforrnaces. For complicated turbulent flows, to 
correctly predict the turbulent structures, a sophisticated non-equilibrium turbulence 
model is needed. 
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